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Abstract 

We present a suite of programs to determine the ground state of the time-independent 
Gross-Pitaevskii equation, used in the simulation of Bose-Einstein condensates. The 
calculation is based on the Optimal Damping Algorithm, ensuring a fast convergence 
to the true ground state. Versions are given for the one-, two-, and three-dimensional 
equation, using either a spectral method, well suited for harmonic trapping poten- 
tials, or a spatial grid. 
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External routines: External FFT or eigenvector routines may be required 

Nature of problem: 

The order parameter (or wave function) of a Bose-Einstein condensate (BEC) is ob- 
tained, in a mean field approximation, by the Gross-Pitaevskii equation (GPE) [1]. 
The GPE is a nonlinear Schrodinger-like equation, including here a confining po- 
tential. The stationary state of a BEC is obtained by finding the ground state of 
the time-independent GPE, i.e., the order parameter that minimizes the energy. In 
addition to the standard three-dimensional GPE, tight traps can lead to effective 
two- or even one-dimensional BECs, so the 2D and ID GPEs are also considered. 
Solution method: 

The ground state of the time-independent of the GPE is calculated using the Opti- 
mal Damping Algorithm [2] . Two sets of programs are given, using either a spectral 
representation of the order parameter [3], suitable for a (quasi) harmonic trapping 
potential, or by discretizing the order parameter on a spatial grid. 
Running time: 

Prom seconds in ID to a few hours for large 3D grids. 
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LONG WRITE-UP 



1 Introduction 



Advances in cooling methods for dilute atomic gases have made it possible to 
attain a new state of matter, the Bose-Einstein condensate (BEC) [Ill2] . As the 
temperature of atoms gets very low, their de Broglie wavelength, an inherently 
quantum character, can become greater than the interatomic distance. At that 
point, bosonic atoms will "condense" into a unique quantum state and become 
indistinguishable parts a of macroscopic quantum object, the BEC. It has now 
been achieved for all stable alkali atoms [5]H1I5]IS|7] . as well as with hydrogen [8j, 
metastable helium PITU] . and for diatomic molecules [TT] . 



Starting from the many-body Hamiltonian describing the cold atoms, it is pos- 
sible to reduce the problem, by considering the order parameter, or wave func- 
tion, for the condensed fraction only. It is governed by a nonlinear Schrodinger 
equation, the Gross-Pitaevskii equation (GPE) 



2m 



Vtrap(x) + A, 



3D 



^(x) = ^V^(x) 



with the normalization condition 



L2 



1, where h is the reduced Planck 



constant, m the mass of the boson, Vtrap a trapping potential spatially confin- 
ing the condensate, and /i the chemical potential of the condensate. Physically, 
the nonlinearity corresponds to the mean field exerted on one boson by all the 
others and is given, for a condensate of bosons in 3D, by 



A 



3D 



fl'3DA^ 



m 



(2) 



The value of a, the scattering length, varies according to the species of bosons 
being considered. The energy associated with the wave function is ob- 
tained according to [T2]|13|14|15lll6j 



N 



|^|V^(x)|2 + Krap( 

2m 



|2 _j_ --^3D 



dx. (3) 



We present here a suite of programs designed to calculate the ground state 
of the GPE, i.e., the order parameter ip{x.) with to the lowest energy. This 
corresponds to the actual condensate order parameter, in the absence of any 
excitation. The problem is thus to find the ground state of the condensate, that 
is a normalized function ipQsi'^) that minimizes E[iIj]. Recall that if Krap is 
continuous and goes to +oo at infinity, and if Asd > 0, the ground state of £^ 
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exists and is unique up to a global phase. In addition, the global phase can be 
chosen such that ipGS is real- valued, and positive on M^. The ground state ipcs 
can be computed using the Optimal Damping Algorithm (ODA), originally 
developed for solving the Hartree-Fock equations [TTPS] . This algorithm is 
garanteed to converge to the ground state. Two different discretizations of the 
order parameter are available in our sets of programs. In one case, a basis 
set of eigenf unctions of the harmonic oscillator is used, which is particularly 
suited for a harmonic (or quasi-harmonic) trapping potential Krap- In this case, 
an efficient method to convert from the spectral representation to a spatial 
grid [in] is employed to treat the nonlinearity. In the other spatial grid 

is used throughout, with the kinetic energy derivative evaluated with the help 
of Fast Fourier Transforms. Note that, in all cases, the value of the energy 
given on output is actually the energy per particle, E[ip]/N. 



2 Optimal Damping Algorithm 

To describe the ODA p!7ffT8] in the context of the GPE, we start by defining 
the operators 

^o = -|^V^ + nx), (4) 
corresponding to the linear part of the GPE ([T]), and 

H{p) =Ho + A3Dp(x) (5) 

the full, nonlinear Hamiltonian, where we have introduced p = {N p[x.) 
is the density of the condensate at point x). 

The ODA is based on the fact that the ground state density matrix 7gs = 
li^Gs) {'4'Gs\ is the unique minimizer of 

inf {£[7], 7 e S{L\R^)), < 7 < /, tr(7) = l} . (6) 

In the above minimization problem, iS(L^(]R^)) denotes the vector space of 
bounded self-adjoint operators on L^(]R^) and / the identity operator on L^(R'^). 
The energy functional £[7] is defined by 

where ^-^(x) = 7(x, x) (7(x, y) being the kernel of the trace-class operator 7). 
The ODA implicitly generates a minimizing sequence 7^ for ([6]), starting, for 
instance, from the initial guess 70 = \4'o){4'o\, where ipo is the ground state 
of Hq. The iterate •jk+i is constructed from the previous iterate 7^ in two steps: 
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Step 1: compute a normalized order parameter ■^^ which minimizes 
f d 

Sk = mf < —£ [(1 - t)-fk + t 



L2 

t=o 



It is easy to check that ip'j^ is in fact the ground state of H{p^^) and that 
either ip'i^ = ^qs (up to a global phase) or < 0. 
• Step 2: compute 

ak = arginf {£: [(1 - t^^ + , ^ e [0, 1]} 

and set 7^+1 = (1 — ak)"fk + Oik\'ipk){'^k\- Note that a can be computed 
analytically, for the function t ^ £ [{1 — t)jk + ^ second order 

polynomial of the form £[yk\ + ^Sfe + yCfe. 

The set 

C = {7 e S{L'{R')), < 7 < tr(7) = 1} 

being convex, 7^ G C for all /c and either 7^ = 7gs or £[7^+1] < S^jk]- In addi- 
tion, it can be proved that, up to a global phase, ip'f. converges to ipcs when k 
goes to infinity. Likewise, pk = p-y,. converges to pes = '4'gs- important to 
note that the sequences ip'j^ and pk can be generated without explicitely com- 
puting 7fc. This is crucial to reduce the overall memory requirement of ODA. 

Let us now describe a practical implementation of ODA, in which only order 
parameters and densities are stored in memory. The algorithm is initialized by 
V'o, from which we derive po = \iI)q\^, /q = (-00, -^oV'o), and ho = {ipo, H{po)ipo). 
The iterations go as follows: 



(1) Calculate the ground state ■^it of H{pk), and p'^ = 

(2) Compute 

fk-i^lHo^'kl 
K^ii^lHipkm, 

K-ii^'k.Hip'kWk)- 

(3) Calculate 

Sk = h'^. — hk, 

Ck = hk + hl-2h', + f',-fk. 

(4) Set ctfe = 1 if Cfc < —Sk, Oik — —Sk/ck otherwise, and 



|2 
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Pk+1 = (1 - ak)Pk + akp'k, 
fk+i = (1 - "fc)/fc + ttfc/fc, 
hk+i = 2-E'opt — fk+i- 

(5) If |sjfc/-Eopt| > £oDA (convergence criterion), go to (1), otherwise compute 
the ground state of H{pk+i), which is the solution sought, and terminate. 

To calculate the ground state of the operators Hq and H{p), the inverse power 
method is used, with the convergence criterion — Ei\ < eip, where E are 
the lowest eigenvalues at consecutive iterations. The inverse power algorithm 
itself uses the conjugated gradient method to solve Hv = u, with u given and 
V unknown. The convergence of the conjugated gradient is controlled by the 
criterion Ecg- The only exception to this is in gpodalDs, where the ground 
states of the operators are found by a matrix eigenproblem solver routine (see 
Sec. OTT]) . 



3 Representations of the GPE 

The Gross-Pitaevskii equation was defined in Eq. ([T]), with the nonlinearity 
Eq. ([2]) in 3D. In this work, we are also considering cases where the con- 
finement Vtrap is so tight in some spatial dimension that the condensate can 
actually be considered as a two-, or even one-dimensional object. This leads 
to different representations of the nonlinearity A and the expression for the 
coupling parameters (720 and Qid can be found in Refs. [20)1211122] . We refer 
to chapter 17 of |2j for a detailed discussion of the validity of the mean field 
approximation in these cases. 



3.1 Spatial grid approach 

If the order parameter is represented on a discretized spatial grid, the calcula- 
tion of the potential energy and the nonlinearity are trivial, as they both act 
locally, while the kinetic energy operator is non-local. By means of a Fourier 
transform, it is possible to convert from position to momentum space, where 
the kinetic operator is local. This is implemented by means of a Fast Fourier 
Transforms (FFTs), allowing to convert back and forth between the two rep- 
resentations, to evaluate each part of the Hamiltonian in the space where it is 
local. 
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3.2 Spectral method 

For many situations, the trapping potential is harmonic, or a close variation 
thereof, z.e.. 



where u is the trapping frequencies in each direction and Vq accounts for 
eventual corrections to a purely harmonic trap. In this case, it is advantageous 
to use a basis set made up of the eigenf unctions of the quantum harmonic 
oscillator. 

We start by rescaling Eq. ([T]), introducing dimensionless lengths {x,y,z), 




ujIx^ + uj^y^ + ujIz^^ + Vq{x, y, z), 



(7) 








and a new order parameter ip defined as 




(9) 



Considering the normalization condition 




(10) 



we take 





such that 



JR3 

The Gross-Pitaevskii equation now reads 






with 




1 



(14) 
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AsD = -TTT^ ^ ^3dA^ = 47raA^ — ^ , (15) 



and 



A^T^- (16) 



Similarly, 



iSW - (17) 
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Using the Galerkin approximation, we can express the order parameter ip as 
a linear combination of a finite number of (orthonormal) basis functions 0, 

Ni Ns 

^(^, y,z) =J2Y.J2 Cijfc0i(x)0j(y)0fc(5), (18) 
i=0 j=0 k=0 

where the (p are chosen as the eigenfunctions of the ID harmonic oscillator, 
i.e., 

Bi|4)*"«)=("4)««- 

In the spectral representation of Eq. ( |T8i) . Eq. ( |T3l) becomes a series of coupled 
equations for the coefficients Qjk, and the first part of the Hamiltonian can be 
evaluated by a simple multiplication, according to Eq. (fT9i) . The second part 
of the Hamiltonian, consisting of the Vq and the nonlinear terms, is local in 
(x, y, z) and couples the different coefficients. Its operation can be calculated 
in a manner similar to what is used for the spatial grid (see Sec. 13.11) : starting 
from the coefficients Cijk, the order parameter ip is evaluated at selected grid 
points (x, ^, S), the local terms are then trivially calculated, and the order 
parameter is transformed back to the spectral representation. This procedure 
can be performed efficiently and accurately using the method described in 

Ref. cni. 



For the 2D case, i.e., when the motion along y is suppressed, we rescale the 
lengths according to Eq. ([8]), which results in 



m\'^/'^ , ,1/4 



A = (^-j {uj.iu^r''' (20) 
for the scaling factor of the order parameter. We thus obtain the 2D GPE 



where 



(21) 

A2D = A2D^ y— I • (22 
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Similarly, we get for the one- dimensional case (where the motion along x and 
y is frozen) 



A 



h 



;;2 



ID 



and 



Aid = Aid 



m 



1/2 



(23) 
(24) 

(25) 



4 Description of the programs 



4-1 gpodaSDg 



This program solves the full 3D GPE ([T]) on a grid. Atomic units are used 
throughout. 



4.1.1 User- supplied routines 



The double precision function potentialV(x,y,z) takes as input the three 
double precision arguments x, y, and z, corresponding to the spatial coordi- 
nates {x,y,z), and returns Vtrs,p{x , y , z) . 

A 3D FFT routine must also be supplied. The program is set up to work with 
the DFFTPACK [23] transform of a real function, and can be linked directly to 
this library. 

If the user wishes to use another FFT, the file f ourierSD . f 90 must be mod- 
ified accordingly. The program first calls fft_init(n), where n is a one- 
dimensional integer array of length 4, the last three elements containing the 
number of grid points in x, y, and z, with the first element corresponding 
to the maximum number of grid points in any direction, i.e., for n(0:3), 
n(0) = maxval (n(l : 3) ) . The program will then call repeatedly the subrou- 
tine f ourier3D(n, fin, f out , direction), with fin and fout double preci- 
sion arrays of dimension (n(l) ,n(2) ,n(3)), and direction an integer. The 
routine should return in array fout the forward Fourier transform of fin if 
direction = 1, and the inverse transform for direction = —1. Any variable 
initialized by f ft_init must be passed to f ourier3D through a module. Note 
that the main program expects to receive the Fourier coefficients (following 
the forward transform) according to: 
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N 

Cl = fn, 
n=l 

m = 2, . . . , 7V/2 + 1 
, m = 2,...,N/2 

where the coefficients correspond to variable f out and the sequence /„ to 
fin. 



4- 1-2 Input parameters 

The input parameters are read from a namelist contained in a file named 
paramsSDg . in, with the following format (the variable type is indicated in 
parenthesis, where dp stands for double precision): 

&params3Dg 

mass = mass of the boson (dp), 

lambda = nonlinearity X^b (dp), 

ng_x = number of grid points in x, (integer) , 

ng_y = number of grid points in y, (integer) , 

ng_z = number of grid points in z, (integer) , 

xmin = first point of the grid in x (dp) , 

xmax = last point of the grid in x (dp) , 

ymin = first point of the grid in y (dp) , 

ymax = last point of the grid in y (dp) , 

zmin = first point of the grid in z (dp) , 

zmax = last point of the grid in z (dp) , 

critODA = convergence criterion for the ODA, Sqda (dp), 

critIP = convergence criterion for the inverse power, ejp (dp) , 

critCG = convergence criterion for the conjugated gradient, scg (dp), 

itMax = maximum number of iterations of the ODA (integer) , 

guess_f rom_f ile = read initial guess from file guessSDg. data? (logical) 
&end 

If the value of the input parameter guess_f rom_f ile is .true., a file named 
guessSDg.data must be present in the local directory. It contains the initial 
guess for the order parameter, and must consist in ng_x x ng_y x ng_z lines, 
each containing the values of the coordinates x, and followed by -0(0;, |/, z). 
Note that the program does not check if the coordinates correspond to the 
grid defined by the input parameters. The program will simply assign the first 
value of %Ij to the first grid point, (xmm, l/min, ^min), then the second value 
to the second grid point in x, with y = T/min and z — Zmin, etc. After rix 



N 



C2 



\m-2 



fn 



COS 



27r(m — l)(n 
TV 



C2m-1 



-Y^fn sin 



n=\ 



27r(m- l)(n- 1) 
N 
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points have been read, the next value of ijj is assigned to the second grid 
point in y, with x = Xmin and z = Zmin, and so on. In other words, the fourth 
column of guessSDg . data contains ip{x, y, z) in standard Fortran format, with 
X corresponding to the first index, y to the second, and z to the third. 

4.1.3 Output files 

The order parameter is written out in file gsSDg.data, with each line contain- 
ing the coordinates x, y, and z, followed by ip{x, y, z). If the algorithm has not 
converged, the file will contain the function obtained at the last iteration. The 
format of gsSDg.data is the same as that of guessSDg. data (see Sec. I4.1.2p . 
such that gsSDg.data can be used as an initial guess for a new run, with for 
instance a different value of A (if the grid is changed, the function must be 
interpolated to the new grid beforehand). 



4.2 gpoda2Dg 

This program solves the 2D GPE on a grid, corresponding to the 3D case 
where motion along y is frozen. Atomic units are used throughout. 

4.2.1 User-supplied routines 

The double precision function potentialV(x,z) takes as input the two double 
precision arguments x and z, corresponding to the spatial coordinates {x,z), 
and returns Vtra.p{x,z). 

A 2D FFT routine must also be supplied. The program is set up to work with 
the DFFTPACK [23] transform of a real function, and can be linked directly to 
this library. For use of another FFT routine, please see Sec. I4.1.1[ 

4.2.2 Input parameters 

The input parameters are read from a namelist contained in a file named 
params2Dg. in. The namelist &params2Dg follows the same format as the 
namelist feparamsSDg presented in Sec. I4.1.2[ with the omission of variables 
ng_y, ymin, and ymax. Also, the parameter lambda corresponds here to g2DN 



If the value of the input parameter guess_f roni_f ile is .true., a file named 
guess2Dg . data must be present in the local directory. The format of the file is 
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similar to that of guess 3Dg.dat a, presented in Sec. 14.1.21 with the exception 
of data corresponding to coordinate y. 

4.2.3 Output files 

The order parameter is written out in file gs2Dg.data, with each line con- 
taining the coordinates x and z, followed by ip{x, z). If the algorithm has not 
converged, the file will contain the function obtained at the last iteration. The 
format of gs2Dg.data is the same as that of guess2Dg.data (see Sec. I4.2.2p . 
such that gs2Dg.data can be used as an initial guess for a new run, with for 
instance a different value of A2D (if the grid is changed, the function must be 
interpolated to the new grid beforehand). 

4.3 gpodalDg 

This program solves the ID GPE on a grid, corresponding to the 3D case 
where motion along x and y is frozen. Atomic units are used throughout. 

4.3.1 User-supplied routines 

The double precision function potentialV(z) takes as input the double pre- 
cision argument z, corresponding to the spatial coordinate z, and returns 

An FFT routine must also be supplied. The program is set up to work with 
the DFFTPACK [23J transform of a real function, and can be linked directly to 
this library. For use of another FFT routine, please see Sec. I4.1.1[ 

4.3.2 Input parameters 

The input parameters are read from a namelist contained in a file named 
paramslDg. in. The namelist feparamslDg follows the same format as the 
namelist &params3Dg presented in Sec. 14.3.21 with the omission of variables 
ng_x, ng_y, xmin, xmax, ymin, and ymax. Also, the parameter lambda corre- 
sponds here to gmN [20j. 

If the value of the input parameter guess_f roni_f ile is .true., a file named 
guess IDg. data must be present in the local directory. It contains the initial 
guess for the order parameter, and must consist in ng_z lines, each containing 
the values of the coordinate z followed by il){z). Note that the program does not 
check if the coordinates correspond to the grid defined by the input parameters. 
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The program will simply assign the first value of ip to the first grid point, -Zmin, 
then the second value to the second grid point in z, and so on. 

4.3.3 Output files 

The order parameter is written out in file gslDg.data, with each line con- 
taining the coordinate z followed by ipi^)- If the algorithm has not converged, 
the file will contain the function obtained at the last iteration. The format of 
gslDg.data is the same as that of guesslDg.data (see Sec. I4.3.2p . such that 
gslDg.data can be used as an initial guess for a new run, with for instance a 
different value of Aid (if the grid is changed, the function must be interpolated 
to the new grid beforehand). 

4.4 gpodaSDs 

This program solves the full 3D GPE (fT3ll using a spectral method. Note that 
the value of mu calculated is actually the rescaled /2 defined by Eq. f|T6l) . 

4-4-^ User-supplied routines 

The double precision function potential VO(x,y,z) takes as input the three 
double precision arguments x, y, and z, corresponding to the rescaled spatial 
coordinates {x,y,z), and returns Vo{x,y,z), defined by Eq. (|T^ . 

4.4-2 Input parameters 

The input parameters are read from a namelist contained in a file named 
paramsSDs . in, with the following format (the variable type is indicated in 
parenthesis, where dp stands for double precision): 

&params3Ds 
lambda = nonlinearity A3D [Eq. ( 173]) / (dp) , 
wxwz = trap frequency ratio uj^/^^z (dp) , 
wywz = trap frequency ratio Uy/uj^ (dp) , 
n_x = highest basis function in x, (integer) , 
n_y = highest basis function in y, Ny (integer) , 
n_z = highest basis function in z, (integer) , 
synimetric_x = symmetric potential in x (logical) , 
syinmetric_y = symmetric potential in y (logical) , 
syinmetric_z = symmetric potential in z (logical) , 
critODA = convergence criterion for the ODA, £oda (dp). 
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critIP = convergence criterion for the inverse power, eip (dp) , 
critCG = convergence criterion for the conjugated gradient, scg (dp) , 
itMax = maximum number of iterations of the ODA (integer), 
guess_f rom_f ile = read initial guess from file guessSDs . data? (logical) 
output_grid = write final order parameter to file gs3Ds_grid.data? (logical) 
&end 

The algorithm used to find the roots of the Hermite polynomial, needed for the 
spectral method pjj, limits the acceptable highest basis function to n < 91. 
The value of the parameters symmetric allow to reduce the size of the basis 
set used, for the case where the additional trapping potential Vq [Eq. ([7])] is 
even along any of the axes. For instance, if Vo{x,y,z) = Vo{—x,y,z), setting 
symmetric_x = .true, will restrict the basis set along x to even functions 
4>{x) [Eq. f|T8|) ]. as the order parameter will present the same parity as the 
trapping potential V^rap- Note that in all cases the parameters n set the index 
of the highest harmonic oscillator eigenfunction used, not the number of basis 
functions used. 

If the value of the input parameter guess_f rom_f ile is .true., a file named 
guessSDs .data must be present in the local directory. It contains the initial 
guess for the order parameter and contains lines with the values of indices 
i, j, and k (all integers), followed by the coefficient Cijk (double precision), 
see Eq. f llSp . If an index is greater than the value of for the corresponding 
spatial axis, or if its parity is not consistent with the chosen symmetry (see 
above), it is ignored. If a set of indices ijk appears more than once, only the 
last value of Cijk is kept, and any Cijk not specified in the file is taken to be 
equal to zero. 

If the value of the input parameter output_grid is .true . , a second namelist 
will be read from the file paramsSDs . in: 

fegridSD 

ng_x = number of grid points in x, (integer) , 
ng_y = number of grid points in y, (integer) , 
ng_z = number of grid points in z, (integer) , 
xmin = first point of the grid in x (dp) , 
xmax = last point of the grid in x (dp) , 
ymin = first point of the grid in y (dp) , 
ymax = last point of the grid in y (dp) , 
zmin = first point of the grid in z (dp) , 
zmax = last point of the grid in z (dp) 
&end 

(see next section for details on usage). 
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4-4-3 Output files 



The order parameter is written out in file gs3Ds .data, with each hne contain- 
ing the indices i, j, and k, followed by the coefficients Cijk of Eq. f[T5]) . If the 
algorithm has not converged, the file will contain the function obtained at the 
last iteration. The format of gs3Ds . data is the same as that of guessSDs . data 
(see Sec. 14.4.21) . such that gsSDs.data can be used as an initial guess for a 
new run, with for instance a different value of A. 

If the value of the input parameter output_grid is . true . , the order param- 
eter is also written out to the file gs3Ds_grid.data, with each line containing 
the coordinates x, y, and z, defined by the namelist &grid3D, followed by 



4-5 gpoda2Ds 

This program solves the a 2D GPE using a spectral method. Note that the 
value of mu calculated is actually the rescaled /i defined by Eq. (|T6l) . 

4-5.1 User-supplied routines 

The double precision function potential VO(x,z) takes as input the three 
double precision arguments x and z, corresponding to the rescaled spatial co- 
ordinates (x, z), and returns Vo{x, z), defined by the 2D equivalent of Eq. f|T4|) . 

4-5.2 Input parameters 

The input parameters are read from a namelist contained in a file named 
params2Ds . in. The namelist &params2Ds follows the same format as the 
namelist &params3Ds presented in Sec. 14.4. 2[ with the omission of variables 
wywz, n_y, and symmetry_y. Also, the parameter lambda corresponds here to 
A2D [Eq. m]- 

If the value of the input parameter guess_f roni_f ile is .true., a file named 
guess2Ds .data must be present in the local directory. The format is the same 
as the file guessSDs .data (Sec. I4.4.2p . except that only indices i and k are 
present. 

If the value of the input parameter output_grid is .true . , a second namelist 
named &grid2D will be read from the file params2Ds . in. This namelist is the 
same as fegridSD of Sec. 14.4.21 without the variables corresponding to y. 
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4-5.3 Output files 

The order parameter is written out in file gs2Ds.data, with a format similar 
to file gs3Ds . data described in Sec. 14.4.3^ except that only indices i and k 
are present. If the value of the input parameter output_grid is .true., the 
order parameter is also written out to the file gs2Ds_grid.data, in the same 
manner as for file gs3Ds_grid . data (Sec. 14.4^31) . but without the y coordinate. 

4-6 gpodalDs 

This program solves the a ID GPE using a spectral method. Note that the 
value of mu calculated is actually the rescaled /2 defined by Eq. (fT6!l . 

4-6.1 User-supplied routines 

The double precision function potential VO(z) takes as input the three double 
precision arguments z, corresponding to the rescaled spatial coordinate z, and 
returns Vo{z), defined by the ID equivalent of Eq. ( fT^ . 

A routine for calculating eigenvalues and eigenvectors must be supplied. The 
program is set up to use the LAPACK |24] routine for the eigenvalue prob- 
lem for a real symmetric matrix. To use another routine, file eigenlD.f90 
has to be modified. The subroutine eigen(n,H,eigenval,eigenvec) takes 
as input the integer n and the double precision array H(n,n). On output, the 
double precision real eigenval and the double precision array eigenvec(n) 
contain repectively the smallest eigenvalue of matrix H and the corresponding 
eigenvector. 

4-6.2 Input parameters 

The input parameters are read from a namelist contained in a file named 
paramslDs . in, with the following format (the variable type is indicated in 
parenthesis, where dp stands for double precision): 

feparamslDs 
lambda = nonlinearity Aid [Eq- ( f^) / (dp) , 
n = highest basis function, N (integer), 
symmetric = spatially symmetric potential? (logical), 
critODA = convergence criterion for the ODA, ^oda (dp), 
itMax = maximum number of iterations of the ODA (integer), 
guess_f rom_f ile = read initial guess from file guesslDs . data? (logical) 
output_grid = write final order parameter to file gslDs_grid.data? (logical) 
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&end 



See Sec. 14.4.21 for restrictions on the value of n and the use of symmetric. 

If the value of the input parameter guess_f rom_f ile is .true., a file named 
guess IDs .data must be present in the local directory. The format is the same 
as the file guesslDs .data (Sec. I4.4.2p . except that only index k is present. 

If the value of the input parameter output_grid is .true . , a second namelist 
named fegridlD will be read from the file paramslDs . in. This namelist is the 
same as fegridlD of Sec. 14.4.21 without the variables corresponding to x and 

y- 

4-6.3 Output files 

The order parameter is written out in file gslDs . data, with a format similar to 
file gslDs.data described in Sec. 14.4. 3[ except that only index k is present. If 
the value of the input parameter output_grid is .true . , the order parameter 
is also written out to the file gslDs_grid.data, in the same manner as for file 
gslDs_grid.data (Sec. 14.4.31) . but without the x and y coordinates. 
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TEST RUN OUTPUT 



Considering a condensate of 10^ ^''Rb atoms, in a harmonic trap of frequency 
<^a; = f-i^y = '^zj \/8 = 27r X 90 Hz with the parameter file paramsSDs . in as 
follows: 

&params3Ds 

lambda = 368. 8d0, 

wxwz = 0.353553390593d0, 

wywz = 0.353553390593d0, 

n_x = 20, 

n_y = 20, 

n_z = 20, 

syinmetric_x = .true., 
symmetric.y = .true., 
syinmetric_z = .true., 
critODA = l.d-8, 
critIP = l.d-8, 
critCG = l.d-8, 
itMax = 100, 

guess_f rom_f lie = .false., 
output _grid = .false. 
&end 

the output will look like: 

GP0DA3DS 

Parameters : 

omega_x / omega_z = . 35355339E+00 
omega_y / omega_z = . 35355339E+00 
Nonlinearity = . 36880000E+03 

Number of basis functions: 11 x 11 x 11 = 1331 
Number of grid points: 41 x 41 x 41 = 68921 

Symmetric in x y z 

Initialization 

Compute the ground state of H_0 
Inverse Power converged in 2 iterations 

— > mu = 0.853553390593000 
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Iteration 1 



Compute the ground state of H(psi_in) 
Inverse Power converged in 27 iterations 

— > mu = 2.19942774785621 
Optimal damping 

slope -22.0705785783271 

step 0.768246751736393 

Eopt 4.08395470599574 



Iteration 66 

Compute the ground state of H(psi_in) 
Inverse Power converged in 2 iterations 

— > mu = 3.90057925938285 
Optimal damping 

slope -1.667024296381214E-008 

step 3.537833144766566E-002 

Eopt 2.87515659549269 



Convergence achieved in 66 iterations 
— > mu = 3.90057925938285 
— > E = 2.87515659549269 



Checking self-consistency 

Inverse Power converged in 2 iterations 

— > mu = 3.90057337542902 
11 norm of psi2out-psi2in: 0.171325E-01 for 68921 grid points 
(0.248582E-06 per grid point) 

(Running time: 14 min on a 2.5 GHz PowerPC G5 Quad.) 
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